perm filename DHSYIQ.SAI[PIC,HE] blob sn#430334 filedate 1979-04-03 generic text, type C, neo UTF8
COMMENT ⊗   VALID 00012 PAGES
C REC  PAGE   DESCRIPTION
C00001 00001
C00003 00002	ENTRY "dhsyiq"
C00008 00003	! This is a short description of what the algorithms really do.
C00012 00004	! this page contains all declarations and global procedures
C00015 00005	! this page contains the macros necessary for the hue calculations
C00018 00006	! this page contains the code to initialize files and related variables
C00021 00007	! this page contains the code to initialize various constants
C00022 00008	! this page contains the code to preload various arrays
C00024 00009	! this page contains the code for beginning the main outer (rows) and inner (columns) loops
C00026 00010	! this page contains the guts of the program
C00029 00011	! this page clears buffers and releases them
C00030 00012	REQUIRE UNSTACK!DELIMITERS
C00031 ENDMK
C⊗;
ENTRY "dhsyiq";
BEGIN "dhsyiq"
REQUIRE "35A" COMPILER!SWITCHES;
REQUIRE "vislib.dcl" SOURCE!FILE;
REQUIRE "⊂⊃<>" DELIMITERS;

!   Rewritten by John R. Kender, July, 1975.
    Modified  by John R. Kender, November, 1976.

This program is a combination and modification of some previous
programs whose purposes were:

1) Transform the original red, green, blue information data into the
psychologically descriptive attributes of intensity, hue, and
saturation.

2) Transform the original red, green, blue information data into the
linear transformations that the color television industry uses.
These are called "Y", "I", and "Q", and correspond roughly to
luminance, orange-red filtration, and magenta filtration,
respectively.

The procedure assumes that the buffer routines have been initialized
by a "bufinit" call outside itself.

Intensity (density) is simply arithmetic mean of the tristimulus
values, rounded.

Hue, saturation and the Y, I, and Q transformations are calculated as
described in the CMU Technical Report, "Saturation, Hue, and Norm-
alized Color:  Calculation, Digitization Effects, and Use" by John R.
Kender.

To save time, many of the repeated calculations are done before-
hand, and their values stored in arrays.  Thus the intensity, satur-
ation, Y, I, and Q values are calculated using table lookups rather
than by straight multiplication or division.  Similarly for part of
the hue calculation.

For ease of programming (mine), the hue calculation is always
performmed assuming the case of a 360 degree output.  If a smaller
byte size is necessary, the values are divided down by 3, 5, or 15.
Note to reader:  feel free to change this.

Restrictions:  The program is designed to handle input files
containing up to and including eight bits per byte.  If it is
necessary to process files with a larger byte size, change the macro
"!biggest!byte!handled!" to the new byte size and recompile.  Some
reprogramming may be necessary for hue at these byte sizes.  Note
also that most of these transformations make little sense if the
input byte size is very small.

Future note:  The performance can be greatly improved if some sort of
anticipatory buffering is included (or else, buffering in units of
blocks).  Note also that the standard picbuf practice of creating a
file of zeros, then reading in a row, modifying it and rewriting it
is unnecessary.  Since the other six files are being created, there
is no need to "update" the new all zero files, and they can be
written directly (say, with a SAIL "arryout" statement).

;
! This is a short description of what the algorithms really do.

Intensity:

    i ← (r + g + b) / 3.0 + .5,  that is, rounded mean.
    As r, g, b integers, integer division with proper prefudge gives
    identical result:
    i ← (sum + 1) % 3, which is easily converted to table lookup.

Hue:

    After deciding which of the six parts of the color triangle to use,
    the arctangent is derived as follows:
    a ← max - min
    b ← mid - min
    That is, a and b are saturated color components.
    ratio ← (a - b) / (a + b)
    This gives relative distance from midpoint of triangle side:
	that is, if ratio=0, angle is 0. ratio=1, angle is 60 degrees.
    Then, ratio ← sqrt(3) * ratio
	gives a proper tangent for the angle under consideration.
    So, angle = arctan (sqrt (3) * (max - mid) / (max - min + mid - min)).
    The value of the tangent is used as an index into a predictor array:
    probe ← atanpredictor[ratio * xfactor]
	where xfactor is determined by the size of the predictor array.
    tantable[probe] is the tangent of the predicted (integer) angle.
        If it is too small, then probe + 1 is used.
    Note that xfactor is a constant, hence ratio * xfactor is really
	(a - b) * (xfactor / (a + b)).  The factor xfactor / (a + b)
	can be converted to table lookup.
    Note also that the multiplication by the coefficient sqrt(3) can
	also be embedded into the tables (as is a factor of 1024
	in the case of tantable).

Saturation:

    s ← (1 - 3 * (r min b min g) / (r + b + g) )
    This gives a saturation in the range 0 to 1.  Then,
    s ← expander * s + .5
    That is, it is plumped up to an integer, rounded.
    Rearranging terms:
    s ← (expander + .5) - ( (3 * expander) / sum) * (r min b min g)
    The term (3 * expander) / sum can be converted to table lookup.

Y, I, and Q:

    The matrix multiplication is:

    y    .509  1.000   .194   r
    i = 1.000 - .460  -.540   g
    q    .403 -1.000   .597   b

    As r, g, b integers, each individual multiplication can be
	accomplished by a table look up.
    Rounding and the proper offsetting of negative values can
	be piggybacked into the array lookup, also.;
! this page contains all declarations and global procedures;

INTERNAL PROCEDURE dhsyiq (STRING ARRAY filenames;
		  BOOLEAN randomize);
BEGIN "dhsyiq!code"
! reassign the following, and recompile, if files are > 8 bits/byte;

DEFINE !biggest!byte!handled! ← ⊂ 8 ⊃ ;

SHORT BOOLEAN continue, no, yes;
SHORT BOOLEAN SAFE ARRAY want[4 : 9];
SHORT INTEGER bytz, colz, f, flag, ii, isrt, jj, jsrt, rowz;
SHORT INTEGER SAFE ARRAY buf[1 : 9], outbytsz[4 : 9];
SHORT INTEGER rp, gp, bp, dp, hp, sp, yp, ip, qp,
              r , g , b , d , h , s , y , i , q ;
SHORT INTEGER achromatic, black, satfact, satrange, sum;
SHORT INTEGER hue!divider, hue!range, hue!rounder, hue!achromatic;
SHORT INTEGER cyan, magenta, yellow;
SHORT INTEGER bytemax, bytemaxx2, bytemaxx3, probe;
SHORT REAL    angle, deg!per!rad, offset!iq, pi, rad!per!deg,
	      ratio, sattop, sqrt3;
SHORT REAL    rrand, grand, brand, sumrand;
STRING dev;

! declarations used for the various arrays;

DEFINE !biggest!value!handled! ← ⊂ 2 ↑ !biggest!byte!handled! - 1 ⊃ ;

SHORT INTEGER SAFE ARRAY intens       [0 : 3 * !biggest!value!handled!];
SHORT INTEGER SAFE ARRAY atanpredictor[0 : 1024                      ];
SHORT REAL    SAFE ARRAY t1024over    [0 : 2 * !biggest!value!handled!];
SHORT REAL    SAFE ARRAY tantable     [0 : 60                        ];
SHORT REAL    SAFE ARRAY satfactover  [0 : 3 * !biggest!value!handled!];
SHORT REAL    SAFE ARRAY yrfact, ygfact, ybfact,
			 irfact, igfact, ibfact,
			 qrfact, qgfact, qbfact
				      [0 :     !biggest!value!handled!];

SIMPLE SHORT REAL PROCEDURE tan (SHORT REAL x);
    RETURN (SIN (x) / COS (x));

! this page contains the macros necessary for the hue calculations;

DEFINE !arctan! (max, mid, min) ← ⊂ (
    IF tantable[
		(probe ← atanpredictor[
				       ratio ← (max - mid)
					       * t1024over[max - min + mid - min]
				      ]
		)
	       ] ≥ ratio THEN
	probe
    ELSE
	probe + 1
) ⊃ ;

! macro for calculating randomized arctan using tables;
DEFINE !arctan!rand! (max, mid, min) ← ⊂ (
    IF tantable[
		(probe ← atanpredictor[
				       ratio ← (max - mid)
					       * 1024 / (max - min + mid - min)
				      ]
		)
	       ] ≥ ratio THEN
	probe
    ELSE
	probe + 1
) ⊃ ;

! macro for exploiting the sixfold hue symmetry;
DEFINE !hue!muncher! (red, green, blue, function) = ⊂ (
IF green > blue THEN
	IF red > green THEN
	    yellow - function (red, green, blue)
	ELSE
	    IF red > blue THEN
		yellow + function (green, red, blue)
	    ELSE
		cyan - function (green, blue, red)
    ELSE
	IF green > red THEN
	    cyan + function (blue, green, red)
	ELSE
	    IF blue > red THEN
		magenta - function (blue, red, green)
	    ELSE
		IF blue > green THEN
		    IF (h ← magenta + function (red, blue, green)) < 360 THEN
			h
		    ELSE
			0
		ELSE
		    IF red > green THEN
			0
		    ELSE
			achromatic
) ⊃ ;

! macros for dividing the calculated hue down;
DEFINE !divided!hue! = ⊂ (
IF (h ← !hue!muncher! (r, g, b, !arctan!)) = achromatic THEN
    hue!achromatic
ELSE
    ((h + hue!rounder) % hue!divider) mod hue!range
) ⊃ ;

DEFINE !divided!hue!rand! = ⊂ (
IF (h ← !hue!muncher! (rrand, grand, brand, !arctan!rand!)) = achromatic THEN
    hue!achromatic
ELSE
    ((h + hue!rounder) % hue!divider) mod hue!range
) ⊃ ;
! this page contains the code to initialize files and related variables;

! read in red, green, and blue files;
FOR f ← 1 thru 3 DO
    BEGIN
	dev ← getdev (filenames[f], "dat");
	indmp (dev, filenames[f], buf[f] ← fndbuf (0), flag);
    END;

! set up switches for files desired;
yes ← -1;
FOR f ← 4 thru 9 DO
    IF LENGTH (filenames[f]) ≠ 0 THEN
	want[f] ← yes;

! set up other file sizes and get buffers for them;
! note that the red file determines the parameters of the others;
rowz ← rows   (buf[1]);
colz ← colms  (buf[1]);
bytz ← bytsz  (buf[1]);
isrt ← isubst (buf[1]);
jsrt ← jsubst (buf[1]);

! are we properly equipped and outfitted;
no ← 0;
IF bytz > !biggest!byte!handled! THEN
    BEGIN
	bprmpt ("You have bitten off more than I can chew.  See comment in source.  Continue? ", continue ← no);
	IF NOT continue THEN RETURN;
    END;

outbytsz[4] ← bytz;
outbytsz[6] ← bytz - 1;
outbytsz[7] ← bytz + 1;
outbytsz[8] ← bytz + 1;
outbytsz[9] ← bytz + 1;

! find the hue byte size and other related variables;
pi ← 3.1415926536;

FOR hue!divider ← 1, 3, 5, 15 DO
    BEGIN
	hue!range ← 360 % hue!divider;
	hue!rounder ← hue!divider % 2;
	outbytsz[5] ← LOG (hue!range) / LOG (2) + 1;
	IF 2 ↑ bytz * 2 * pi ≥ hue!range THEN
	    DONE;
    END;

hue!achromatic  ← hue!range + 1;  ! note: convention only;

! initialize output files;
FOR f ← 4 thru 9 DO
    IF want[f] THEN
	BEGIN
	    getbuf (rowz, colz, outbytsz[f], buf[f] ← fndbuf (0));
	    putsub (isrt, jsrt, buf[f]);
	END;
! this page contains the code to initialize various constants;

black       ← 0;    ! note: convention only (and a bad one at that);
bytemax     ← 2 ↑ bytz - 1;
bytemaxx2   ← bytemax * 2;
bytemaxx3   ← bytemax * 3;
deg!per!rad ← 180 / pi;
offset!iq   ← 2 ↑ (outbytsz[8] - 1);
rad!per!deg ← pi / 180;
sqrt3       ← sqrt (3);
satrange    ← 2 ↑ outbytsz[6] - 1;
satfact     ← satrange * 3;
sattop      ← satrange + .5;

! code for special hue stuff;
cyan    ← 180;
magenta ← 300;
yellow  ← 60;
achromatic ← 361;
! this page contains the code to preload various arrays;

! intensity;
FOR ii ← 0 thru bytemaxx3 DO
    intens[ii] ← (ii + 1) % 3;  ! includes rounding;

! hue;
FOR ii ← 0 thru 60 DO
    tantable[ii] ← 1024 * (tan ((ii + .5) * rad!per!deg) / sqrt3);

FOR ii ← 0 thru 1024 DO
    atanpredictor[ii] ← deg!per!rad * atan (ii * sqrt3 / 1024) + .5;

t1024over[0] ← 0;  ! never accessed;
FOR ii ← 1 thru bytemaxx2 DO
    t1024over[ii] ← 1024 / ii;

! saturation;
satfactover[0] ← 0;  ! never accessed, as sum = 0 is the singularity;
satfactover[1] ← 0;  ! if sum = 1, then (r min g min b) = 0;
satfactover[2] ← 0;  ! if sum = 2, then (r min g min b) = 0;
FOR ii ← 3 thru bytemaxx3 DO
    satfactover[ii] ← satfact / ii;

! y, i, and q;
FOR ii ← 0 thru bytemax DO
    BEGIN
	yrfact[ii] ←   .509 * ii                 ;
	ygfact[ii] ←  1.000 * ii                 ;
	ybfact[ii] ←   .194 * ii + .5            ; ! includes rounding;
	irfact[ii] ←  1.000 * ii                 ;
	igfact[ii] ← - .460 * ii                 ;
	ibfact[ii] ← - .540 * ii + .5 + offset!iq; ! includes rounding;
	qrfact[ii] ←   .403 * ii                 ;
	qgfact[ii] ← -1.000 * ii                 ;
	qbfact[ii] ←   .597 * ii + .5 + offset!iq; ! includes rounding;
    END;
! this page contains the code for beginning the main outer (rows) and inner (columns) loops;

FOR ii ← 1 thru rowz DO
    BEGIN "ii"
	rp ← inptr  (ii, 1, buf[1]);
	gp ← inptr  (ii, 1, buf[2]);
 	bp ← inptr  (ii, 1, buf[3]);
	IF want[4] THEN
	    dp ← outptr (ii, 1, buf[4]);
	IF want[5] THEN
	    hp ← outptr (ii, 1, buf[5]);
	IF want[6] THEN
	    sp ← outptr (ii, 1, buf[6]);
	IF want[7] THEN
	    yp ← outptr (ii, 1, buf[7]);
	IF want[8] THEN
	    ip ← outptr (ii, 1, buf[8]);
	IF want[9] THEN
	    qp ← outptr (ii, 1, buf[9]);
	FOR jj ← 1 thru colz DO
	    BEGIN "jj"
		r ← ILDB (rp);
		g ← ILDB (gp);
		b ← ILDB (bp);
		IF want[4] OR (want[6] AND NOT randomize) THEN
		    sum ← r + g + b;
		IF randomize THEN
		    BEGIN
			rrand ← r + RAN (0);
			grand ← g + RAN (0);
			brand ← b + RAN (0);
			IF want[6] THEN
			    sumrand ← rrand + grand + brand;
		    END;
! this page contains the guts of the program;

		! intensity;
		IF want[4] THEN 
		    IDPB (intens[sum] 
			, dp);
		! hue;
		IF want[5] THEN 
		    IF NOT randomize THEN
			IF hue!divider = 1 THEN
			    IDPB (!hue!muncher! (r, g, b, !arctan!)
				, hp)
			ELSE
			    IDPB (!divided!hue!
				, hp)
		    ELSE
			IF hue!divider = 1 THEN
			    IDPB (!hue!muncher! (rrand, grand, brand, !arctan!rand!)
				, hp)
			ELSE
			    IDPB (!divided!hue!rand!
				, hp);
		! saturation;
		IF want[6] THEN
		    IF NOT randomize THEN
			IDPB (IF sum = 0 THEN
				black
			    ELSE
				sattop -  (r MIN g MIN b) * satfactover[sum]
				! note that the type of the above expression is
				  integer, due to the "THEN" part;
			    , sp)
		    ELSE
			IDPB (IF sumrand = 0 THEN
				black
			    ELSE
				satrange * (1 - 3 * (rrand MIN grand MIN brand) / sumrand) + .5
			    , sp);
		! y;
		IF want[7] THEN 
		    IDPB (y ← yrfact[r] + ygfact[g] + ybfact[b]
			, yp);

		! i;
		IF want[8] THEN 
		    IDPB (i ← irfact[r] + igfact[g] + ibfact[b]
			, ip);

		! q;
		IF want[9] THEN 
		    IDPB (q ← qrfact[r] + qgfact[g] + qbfact[b]
			, qp);

	    END "jj";

	IF (ii MOD 25) = 0 THEN
	    outst (CVS (ii) & " Rows happily munched so far!" & crlf);
    END "ii";
! this page clears buffers and releases them;

FOR f ← 1 thru 3 DO
    frebuf (buf[f]);
FOR f ← 4 thru 9 DO
    IF want[f] THEN
	BEGIN
	    dev ← getdev (filenames[f], "dat");
	    outdmp (dev, filenames[f], buf[f], flag);
	    frebuf (buf[f]);
	END;

END "dhsyiq!code";
REQUIRE UNSTACK!DELIMITERS;
END "dhsyiq";